reading the functions in that i’ll use

source(here::here('run_standard_deseq.R'))
source(here::here('make_volcano_plot.R'))

reading in the files

feature_counts_file_path <- file.path(here::here("data","feature_counts_bdnf_No4"))

running a deseq 1hr

one_hour_dds <- run_standard_deseq(feature_counts_file_path, base_grep = "CONTROL",
                              contrast_grep = "BDNF", 
                              grep_pattern = "1hr",
                              baseName = "control",
                              contrastName = 'BDNF')
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## The following object is masked from 'package:Biobase':
## 
##     rowMedians
## 
## Attaching package: 'data.table'
## The following object is masked from 'package:SummarizedExperiment':
## 
##     shift
## The following object is masked from 'package:GenomicRanges':
## 
##     shift
## The following object is masked from 'package:IRanges':
## 
##     shift
## The following objects are masked from 'package:S4Vectors':
## 
##     first, second
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## The following object is masked from 'package:purrr':
## 
##     transpose
## [1] "BDNF_1B_1hr"    "BDNF_2B_1hr"    "BDNF_3B_1hr"    "CONTROL_1A_1hr"
## [5] "CONTROL_2A_1hr" "CONTROL_3A_1hr"
## [1] "BDNF_1B_1hr"    "BDNF_2B_1hr"    "BDNF_3B_1hr"    "CONTROL_1A_1hr"
## [5] "CONTROL_2A_1hr" "CONTROL_3A_1hr"
## [1] "Filtered Genes by CPM greater than 0.5 in a least 2 samples"
## keep
## FALSE  TRUE 
## 35930 24739
## converting counts to integer mode
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## 
## out of 24739 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 297, 1.2%
## LFC < 0 (down)     : 145, 0.59%
## outliers [1]       : 0, 0%
## low counts [2]     : 11511, 47%
## (mean count < 58)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## Joining, by = "Geneid"

make a volcano plot 1hr

one_hour_volcano <- make_volcano_plot(one_hour_dds$results_table)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

deseq 2hr

two_hour_dds <- run_standard_deseq(feature_counts_file_path, base_grep = "CONTROL",
                              contrast_grep = "BDNF", 
                              grep_pattern = "2hr",
                              baseName = "control",
                              contrastName = 'BDNF')
## [1] "BDNF_1D_2hr"    "BDNF_2D_2hr"    "BDNF_3D_2hr"    "CONTROL_1C_2hr"
## [5] "CONTROL_2C_2hr" "CONTROL_3C_2hr"
## [1] "BDNF_1D_2hr"    "BDNF_2D_2hr"    "BDNF_3D_2hr"    "CONTROL_1C_2hr"
## [5] "CONTROL_2C_2hr" "CONTROL_3C_2hr"
## [1] "Filtered Genes by CPM greater than 0.5 in a least 2 samples"
## keep
## FALSE  TRUE 
## 27648 33021
## converting counts to integer mode
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## 
## out of 33021 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 1420, 4.3%
## LFC < 0 (down)     : 6390, 19%
## outliers [1]       : 0, 0%
## low counts [2]     : 14085, 43%
## (mean count < 20)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## Joining, by = "Geneid"

volcano 2hr

two_hr_volcano <- make_volcano_plot(two_hour_dds$results_table)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

deseq 6h

six_hr_dds <- run_standard_deseq(feature_counts_file_path, base_grep = "CONTROL",
                              contrast_grep = "BDNF", 
                              grep_pattern = "6hr",
                              baseName = "control",
                              contrastName = 'BDNF')
## [1] "BDNF_1F_6hr"    "BDNF_2F_6hr"    "BDNF_3F_6hr"    "CONTROL_1E_6hr"
## [5] "CONTROL_2E_6hr" "CONTROL_3E_6hr"
## [1] "BDNF_1F_6hr"    "BDNF_2F_6hr"    "BDNF_3F_6hr"    "CONTROL_1E_6hr"
## [5] "CONTROL_2E_6hr" "CONTROL_3E_6hr"
## [1] "Filtered Genes by CPM greater than 0.5 in a least 2 samples"
## keep
## FALSE  TRUE 
## 27726 32943
## converting counts to integer mode
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## 
## out of 32943 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 963, 2.9%
## LFC < 0 (down)     : 2487, 7.5%
## outliers [1]       : 2, 0.0061%
## low counts [2]     : 17245, 52%
## (mean count < 31)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## Joining, by = "Geneid"

volcano 6h

six_hr_volcano <- make_volcano_plot(six_hr_dds$results_table)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

now we are going to use the deseq data to do gene ontology hour 1

results_1hr <- one_hour_dds$results_table

up_1hr_ensembl <- results_1hr %>% 
  dplyr::filter(padj < 0.05 & log2FoldChange > 0 ) %>% 
  mutate(ensembl = gsub("\\..*", "", Geneid)) %>% 
  arrange(-log2FoldChange) %>% 
  pull(ensembl)

down_1hr_ensembl <- results_1hr %>% 
  dplyr::filter(padj < 0.05 & log2FoldChange < 0 ) %>% 
  mutate(ensembl = gsub("\\..*", "", Geneid)) %>% 
  arrange(log2FoldChange) %>% 
  pull(ensembl)

gp_up_1hr = gost(up_1hr_ensembl, organism = "hsapiens", ordered_query = TRUE) 
gostplot(gp_up_1hr)
gp_down_1hr = gost(down_1hr_ensembl, organism = "hsapiens", ordered_query = TRUE)
gostplot(gp_down_1hr)

GO hr 2

results_2hr <- two_hour_dds$results_table

up_2hr_ensembl <- results_2hr %>% 
  filter(padj < 0.05 & log2FoldChange > 0 ) %>% 
  mutate(ensembl = gsub("\\..*", "", Geneid)) %>% 
  pull(ensembl)

down_2hr_ensembl <- results_2hr %>% 
  filter(padj < 0.05 & log2FoldChange < 0 ) %>% 
  mutate(ensembl = gsub("\\..*", "", Geneid)) %>% 
  pull(ensembl)



gp_up_2hr <-  gost(up_2hr_ensembl, organism = "hsapiens")

gp_down_2hr <-  gost(down_2hr_ensembl, organism = "hsapiens")

go 6h

results_6hr <- six_hr_dds$results_table

up_6hr_ensembl <- results_6hr %>% 
  filter(padj < 0.05 & log2FoldChange > 0 ) %>% 
  mutate(ensembl = gsub("\\..*", "", Geneid)) %>% 
  pull(ensembl)

down_6hr_ensembl <- results_6hr %>% 
  filter(padj < 0.05 & log2FoldChange < 0 ) %>% 
  mutate(ensembl = gsub("\\..*", "", Geneid)) %>% 
  pull(ensembl)

gp_up_6hr = gost(up_6hr_ensembl, organism = "hsapiens")
gp_down_6hr = gost(down_6hr_ensembl, organism = "hsapiens")
bg_ids_1hr <- one_hour_dds$results_table$Geneid %>% gsub("\\..*", "", .) 

e_go_cc_1hr_up <- enrichGO(gene          = up_1hr_ensembl,
                universe      = bg_ids_1hr,
                keyType = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "CC",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

e_go_mf_1hr_up <- enrichGO(gene          = up_1hr_ensembl,
                universe      = bg_ids_1hr,
                keyType = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "MF",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

e_go_bp_1hr_up <- enrichGO(gene          = up_1hr_ensembl,
                universe      = bg_ids_1hr,
                keyType = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

if(any(e_go_cc_1hr_up@result$p.adjust <0.01)){
  cnetplot(e_go_cc_1hr_up)
  dotplot(e_go_cc_1hr_up)
  }

if(any(e_go_bp_1hr_up@result$p.adjust <0.01)){
  cnetplot(e_go_bp_1hr_up)
  dotplot(e_go_bp_1hr_up)}

bg_ids_2hr <- two_hour_dds$results_table$Geneid %>% gsub("\\..*", "", .) 

e_go_cc_2hr_up <- enrichGO(gene          = up_2hr_ensembl,
                universe      = bg_ids_2hr,
                keyType = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "CC",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

e_go_mf_2hr_up <- enrichGO(gene          = up_2hr_ensembl,
                universe      = bg_ids_2hr,
                keyType = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "MF",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

e_go_bp_2hr_up <- enrichGO(gene          = up_2hr_ensembl,
                universe      = bg_ids_2hr,
                keyType = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

if(any(e_go_cc_2hr_up@result$p.adjust <0.01)){
  cnetplot(e_go_cc_2hr_up)
  dotplot(e_go_cc_2hr_up)
  }

if(any(e_go_bp_2hr_up@result$p.adjust <0.01)){
  cnetplot(e_go_bp_2hr_up)
  dotplot(e_go_bp_2hr_up)}

bg_ids_6hr <- six_hr_dds$results_table$Geneid %>% gsub("\\..*", "", .) 

e_go_cc_6hr_up <- enrichGO(gene          = up_6hr_ensembl,
                universe      = bg_ids_6hr,
                keyType = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "CC",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

e_go_mf_6hr_up <- enrichGO(gene          = up_6hr_ensembl,
                universe      = bg_ids_6hr,
                keyType = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "MF",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

e_go_bp_6hr_up <- enrichGO(gene          = up_6hr_ensembl,
                universe      = bg_ids_6hr,
                keyType = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)

if(any(e_go_cc_6hr_up@result$p.adjust <0.01)){
  cnetplot(e_go_cc_6hr_up)
  dotplot(e_go_cc_6hr_up)
  }

if(any(e_go_bp_6hr_up@result$p.adjust <0.01)){
  cnetplot(e_go_bp_6hr_up)
  dotplot(e_go_bp_6hr_up)}